home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsI / Original / SquareRoot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  2.0 KB  |  82 lines  |  [TEXT/MPS ]
  1. /*
  2. A High Speed, Low Precision Square Root
  3. by Paul Lalonde and Robert Dawson
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7.  
  8. /* SPARC implementation of a fast square root by table 
  9.  * lookup.
  10.  * SPARC floating point format is as follows:
  11.  *
  12.  * BIT 31     30     23     22     0
  13.  *     sign    exponent    mantissa
  14. */
  15.  
  16. #include <math.h>
  17. static short sqrttab[0x100];    /* declare table of square roots */
  18. void build_table()
  19. {
  20.     unsigned short i;
  21.     float f;
  22.     unsigned int *fi=&f;    /* to access the bits of a float in  */
  23.                             /* C quickly we must misuse pointers */
  24.     for(i=0; i<= 0x7f; i++)
  25.     {
  26.         *fi = 0;
  27.  
  28.         /* Build a float with the bit pattern i as mantissa
  29.          * and an exponent of 0, stored as 127
  30.           */
  31.  
  32.         *fi = (i << 16) | (127 << 23);
  33.         f = sqrt(f);
  34.  
  35.         /* Take the square root then strip the first 7 bits of
  36.          * the mantissa into the table
  37.         */
  38.  
  39.         sqrttab[i] = (*fi & 0x7fffff) >> 16;
  40.  
  41.         /* Repeat the process, this time with an exponent of
  42.          * 1, stored as 128
  43.             */
  44.  
  45.  
  46.         *fi = 0;
  47.         *fi = (i << 16) | (128 << 23);
  48.         f = sqrt(f);
  49.         sqrttab{i+0x80] = (*fi & 0x7fffff) >> 16;
  50.     }
  51. }
  52.  
  53. /*
  54.  * fsqrt - fast square root by table lookup
  55. */
  56.  
  57. float fsqrt(float n)
  58. {
  59.     unsigned int *num = &n;    /* to access the bits of a float in C
  60.                              * we must misuse pointers */
  61.                             
  62.     short e;                    /* the exponent */
  63.     if (n == 0) return (0);    /* check for square root of 0 */
  64.     e = (*num >> 23) - 127;    /* get the exponent - on a SPARC the */
  65.                             /* exponent is stored with 127 added */
  66.     *num & = 0x7fffff:        /* leave only the mantissa */
  67.     if (e & 0x01) *num | = 0x800000;
  68.                             /* the exponent is odd so we have to */
  69.                             /* look it up in the second half of  */
  70.                             /* the lookup table, so we set the   */
  71.                             /* high bit                   */
  72.     e >>= 1:                    /* divide the exponent by two */
  73.                             /* note that in C the shift */
  74.                             /* operators are sign preserving */
  75.                             /* for signed operands */
  76.     /* Do the table lookup, based on the quaternary mantissa,
  77.      * then reconstruct the result back into a float
  78.      */
  79.     *num = ((sqrttab[*num >> 16]) << 16) | ((e + 127) << 23);
  80.     return(n);
  81. }    
  82.